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Tunneling density of states measurements of disordered superconducting (SC) Al films in high 
Zeeman fields reveal a significant population of subgap states which cannot be explained by stan- 
dard BCS theory. We provide a natural explanation of these excess states in terms of a novel 
disordered Larkin-Ovchinnikov (dLO) phase that occurs near the spin-paramagnetic transition at 
the Chandrasekhar-Clogston critical field. The dLO superconductor is characterized by a pairing 
amplitude that changes sign at domain walls. These domain walls carry magnetization and support 
Andreev bound states that lead to distinct spectral signatures at low energy. 



A central theme in condensed matter physics is the 
quest for new states of matter with unusual arrangements 
of interacting electrons, spins, and atoms. The interplay 
between superconductivity and magnetism is an espe- 
cially rich source of interesting physics that gives rise to 
various types of exotic superconductors such as cupratcs, 
pnictides, ruthenates, and heavy-fermion materials^'^. 
There is also, however, the possibility of exotic super- 
conductivity of a different type, which arises when a con- 
ventional BCS superconductor at low temperature is sub- 
jected to an external Zeeman field. In the simplest sce- 
nario, the superconductor undergoes a first-order tran- 
sition into a polarized normal Fermi liquid^'^ when the 
Zeeman splitting becomes of the order of the supercon- 
ducting gap Ao at the Chandrasekhar-Clogston critical 
field ijlbHcc ~ Ao/\/2. However, nature has a more 
intriguing way of resolving the tussle: the electrons can 
self-organize into a novel intermediate state known as 
a Fuldc-Fcrrcll-Larkin-Ovchinnikov (FFLO) state near 
Hcc-^^^° An FFLO state consists of regions of posi- 
tive and negative pairing amplitude separated by do- 
main walls where the magnetization is piled up; it can 
be thought of as an "electronic liquid crystal," an ex- 
ample of emergent microscale phase separation. Inter- 
est in FFLO physics crosses traditional boundaries be- 
tween condensed matter, cold atomic gases^^, quantum 
chromodynamics^^, nuclear physics, and astrophysics^^, 
and there is currently an intense effort to search for FFLO 
phases in superconductors as well as in cold atoms^^. 

Hitherto, only thermodynamic signatures of the FFLO 
phase have been reported, and these have been limited 
to a few layered organic superconductors and the heavy 
fermion material CeCoIn5^^~^^. The realization of FFLO 
in traditional superconducting systems has been ham- 
pered by its sensitivity to disorder and spin-orbit scat- 
tering. Notwithstanding these issues, we show that even 
in the presence of disorder, where the fully coherent 
FFLO phase is suppressed, spectroscopic manifestations 
of FFLO fluctuations are readily observable. 

Main results: We present density of states (DoS) 
calculations based on a disordered attractive Hubbard 



model, along with low-temperature tunneling DoS mea- 
surements on ultra-thin Al films. We show that, con- 
trary to popular belief, FFLO physics is not completely 
washed out by disorder. In fact, over a significant range 
of Zeeman fields we find a disordered Larkin-Ovchinnikov 
(dLO) state characterized by bound states in domain 
walls and low-energy spectral weight, which provides a 
natural explanation of the experimental anomalies.^* Our 
calculations self-consistently account for the disorder and 
allow the pairing amplitude to adjust to the disorder pro- 
file. The novel dLO phase is robust to variations in field 
and disorder, and imprints a unique signature in the low- 
energy DoS within the superconducting gap. 

Experimental setup: In the present study planar tun- 
nel junctions formed on 3 nm-thick Al films were used 
to extract the low temperature quasiparticlc DoS. Alu- 
minum has a well documented low spin-orbit scatter- 
ing rate^° and superconducting transition temperature 
Tc = 2.7 K with a zero field gap Ao « 0.43 mV in 
thin film form. [For sample preparation see supplement]. 
Measurements of resistance and tunneling were carried 
out on an Oxford dilution refrigerator using a standard 
dc four-probe technique. Magnetic fields of up to 9 T 
were applied using a superconducting solenoid. A me- 
chanical rotator was cimployed to orient the sample in 
situ with a precision of ^ 0.1°. The films were mod- 
erately disordered with sheet resistances of the order of 
Ifcil, well below the quantum of resistance for supercon- 
ductivity Rq = h/Ae^ = 6.4 kO. 

Experimental results and comparison with standard 
BCS theory: We present measurements of the tun- 
neling conductance G of Al films, which is mainly pro- 
portional to the superconducting DoS at the low tem- 
peratures used. Figure 1(a) shows the bias dependence 
G{V) in a parallel field H = 4.75 T at 100 mK, in which 
the BCS coherence peaks have been Zeeman-split by the 
applied field. Figure 1(b) shows the parallel-field de- 
pendence of the zero-bias tunneling conductance G(0), 
which is zero in the conventional superconducting state 
[H < Hq w 2.8 T) and constant in the normal state 
{H > Hc\\ ~ 6.1 T); however, there is a significant tail 
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FIG. 1: (a) Tunneling conductance G{V) normalized by normal state conductance G„ ~ (1 ki7)^^ for a 24 A superconducting 
Al film in a 4.75 T parallel field at 100 mK (symbols=experiment, curve=homogeneous theory), (b) Zero-bias tunneling 
conductance G(0) at 60 mK as a function of parallel field H. Between Hq ~ 2.8 T and H^w ~ 6.1 T, the homogeneous 
theory (blue curve) significantly underestimates the number of states near the Fermi energy, and even when the temperature 
is artificially increased (red curve) it is unable to describe the broad tail in G{0). We ascribe the discrepancy to a disordered 
LO phase. (Inset) Tunnel conductance as a function of H± = 4.5sin(0) where 6 is the tilt angle 9. The solid lines are a linear 
least-squares fit to the data. The sharp V-shaped minimum allows us to accurately determine parallel alignment. 



in G(0) over a range of fields Ho < H < Hc\\- The 
colored curves in Fig. 1(a) and 1(b) are obtained within 
homogenous BCS mean field theory by solving the Usadel 
equations for the disorder-averaged semiclassical Green's 
functions together with the self-consistent equations for 
the uniform order parameter and the internal magnetic 
field. The parameters involved are the gap energy, spin- 
orbit scattering rate, the orbital depairing rate, and the 
antisymmetric Fermi-liquid parameter; they are deter- 
mined by fits^^'^^ to full spectra as in Fig. 1(a). 

The observed excess zero-bias conductance G(0) can 
have various origins, (i) Imperfect alignment: The inset 
of Fig. 1(b), shows G(0) at several alignment angles be- 
tween the film plane and the applied field. It is evident 
that our alignment mechanism is precise enough to find 
parallel orientation within the limits of the sensitivity of 
the tunneling conductance to H± , the perpendicular field 
component, (ii) Junction leakage is ruled out because all 
of the junctions used in this study had a very low zero- 
bias conductance in zero field, G(2mV)/G(0) ~ 10^ — 10^ 
at 100 mK. (iii) Material inhomogeneities: In princi- 
ple could lead to broadened transitions, however, the 
zero-field gap in Al (and hence the nominal critical field 
hcc) varies by only 20% over a very wide range of sheet 
resistance^'^ and averaging over a distribution of gaps fails 
to explain the large range of ify over which G(0) is finite, 
(iv) Pair-breaking: These effects scale as Dd^, where D 
is the normal state diffusivity and d is the film thick- 



ness. For our films as d is decreased from 3 nm to 2 nm, 
D decreases by an order of magnitude, but G(G) hardly 
changes. Furthermore, recent tunneling measurements 
of Al-EuS bilayers have shown that a comparable G(0) is 
produced by an interface-induced exchange field, which 
is a pure Zeeman field with no orbital depairing effects.^'' 

Disordered LO states and excess low- energy spectral 
weight: Having ruled out all the above explanations, 
we now argue that the anomalous excess zero-bias con- 
ductance at intermediate fields is an intrinsic property 
of the condensate due to the development of an exotic 
disordered Larkin-Ovchinnikov (dLO) phase with an in- 
homogeneous pairing amplitude and magnetization. 

Our model consists of the attractive Hubbard Hamil- 
tonian with a disorder potential and a Zeeman field, 

rr'a ra 

-|C/|^(nrt-i)(7M-^) (1) 
r 

where t^r' are hopping amplitudes (equal to t, taken as 
the unit of energy) between nearest-neighbor sites r and 
r', TLra = c\^Cj.^ is the number operator for fermions of 
spin index a = ±1 at site r, ^ is the average chemical 
potential, h is the Zeeman field, and U is the local pair- 
wise Hubbard interaction. The disorder potential Vr at 
each site is picked independently from a uniform distri- 
bution on [— We calculate the local densities n-ra, 
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FIG. 2: Root-mean-square pairing amplitude Arms, average 
magnetization niavg, and Fermi-level density of states A''(0) as 
functions of Zeeman field h, in units of the hopping amplitude 
t (see Eq. (1)). For hd < h < hc2 there is a disordered LO 
state with coexistent pairing and magnetization, in which the 
gap is partially filled in. The results are obtained using BdG 
simulations on a 36 x 36 Hubbard model at weak disorder 
W = It (well below the critical disorder^^ for the destruction 
of superconductivity Wc ~ 3t), nonzero chemical potential 
/i = — 0.25t to avoid perfect nesting effects at half-filling, low 
temperature T — O.lt, and a relatively large attraction \U\ = 
At so that the coherence length is less than the system size. 
h = ^gfiBH, where <; ~ 2 is the (j-factor, /is is the Bohr 
magneton, and H is the parallel field. 

pairing amplitude Ar — \U\ (cr^Cr-f), and spin-dependent 
DoS N„{E) within a fully self-consistent Bogoliubov-de 
Gennes (BdG) framework including all Hartree shifts (see 
supplement for details) . A phase diagram for this system 
was obtained in Ref. 25; in this paper we focus on spectral 
features. 

As illustrated in Fig. 2, if A is restricted to be uni- 
form, we find that the BCS- and normal-state free en- 
ergies cross at hcc = l-Olt, the critical field for the 
first-order Chandrasekhar-Clogston transition (here hcc 
differs from Ao/\/2 due to the moderate value of U). 
However, if A(r) is allowed to be inhomogeneous, BdG 
calculations predict two transitions, at a lower critical 
field hci = 0.85i and an upper critical field hc2 = 1.75i. 
The intermediate state {hd < hcc < ^c2) has both a 
finite pairing amplitude and a finite magnetization. 

A physical understanding is provided in Fig. 3, which 
shows the local pairing amplitude A(r), local magneti- 
zation m{r) — i ~ "-^(i")]! s-iid spatially averaged 
DoS's of up and down spins Na-{E), for various values 
of h. At low fields the system is a BCS superconductor 



with a nearly uniform order parameter A(r) « Aq, whose 
DoS contains coherence peaks at ±A ± h slightly broad- 
ened by inhomogeneous Hartree shifts^^'^^. At high fields 
the system is normal (non-superconducting) with nearly 
uniform magnetization. At intermediate fields there is a 
disordered Larkin-Ovchinnikov (dLO) state with the fol- 
lowing features: (1) There is a strong modulation of the 
pairing amplitude A(r) which changes sign between pos- 
itive and negative values. The oscillations at wavevector 
qlo ~ 2fci? are partially disrupted by the disorder poten- 
tial. (2) The magnetization is finite in the domain walls 
where the pairing amplitude is small. (3) There is signif- 
icant low-energy weight in the DoS, as illustrated in the 
rightmost column of Fig. 3. This is the main new result 
of this paper, and it is a likely explanation for the similar 
low-energy weight seen in experiments (Fig. 1). 

Origin of low-energy states: When the Zeeman field 
exceeds a certain lower critical field, magnetization be- 
gins to penetrate the sample in the form of domain walls 
(brown regions in Fig. 4(a)). The majority electrons 
are unable to enter the superconducting regions due to 
the gap, and so they are confined to the domain walls 
by Andreev reflection, forming Andreev bound states 
with a distribution of energies. Whereas in a clean LO 
state^'^''^^ tunneling between domain walls gives rise to 
subgap bands, in a dLO state the bound states are likely 
to remain localized, but they still contribute to the low- 
energy DoS. Indeed, comparing Figs. 4(d) and (e) shows 
that the low-energy weight is concentrated in the same 
regions as the magnetization. The tunneling DoS (unlike 
transport measurements) is sensitive to local electronic 
structure, and hence the low-energy spectral signatures 
of LO should remain even when phase fluctuations pre- 
vent the development of long-range LO order. 

We conclude that dLO physics is a likely explanation 
of the longstanding mystery of excess zero-bias tunnel- 
ing conductance of Al films near the spin-paramagnctic 
transition.^® Our results suggest that the parallel- ficld- 
tuned^"''^^ superconductor-insulator transition (SIT) oc- 
curs via a dLO phase in which the gap is filled in by 
Andreev bound states. This scenario is distinct from the 
zero-field thickness-tuned "fermionic" SIT where the gap 
c/oses32-34^ and from the "bosonic" SIT26.27,35-37 ^^^^^ 

the gap appears to remain finite across the SIT. 
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FIG. 3: The first two columns show spatial maps of the local pairing amplitude A and the magnetization m. The third 
column show the densities of states (DoS's) of up and down electrons Ncr{E). The last column shows the total DoS N{E). 
For intermediate fields (e.g., h/t = 0.95 and h/t = 1.2) the system exhibits disordered Larkin-Ovchinnikov states with domain 
walls at which m is finite, A changes sign, and the DoS becomes finite at low energy. Other parameters are as in Fig. 2. 
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FIG. 4: (a) Combined plot of m(r) and A(r) for h/t — 1 (other parameters as in Fig. 2). Red (blue) indicates regions where 
A(r) is large and positive (negative). Brown regions, where the magnetization m(r) is large, occur at domain walls where A 
changes sign. White regions are hills or valleys of the disorder potential corresponding to empty sites or localized pairs that 
participate in neither superconductivity nor magnetism, (b) and (c) show oscillations of A along the vertical dashed line in panel 
(a), (d) and (e) show the correspondence between magnetization m(r) and low-energy spectral weight /(r) = /^g^jj dE Nr{E). 
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Mystery of Excess Low Energy States 
in a Disordered Superconductor in a Zeeman Field: 
Supplementairy Information 



I. SAMPLE PREPARATION 

In the present study tunnel junctions were formed by 
first depositing a 3 nm thick Al film e-beam deposition 
of 99.999% Al stock onto fire polished glass microscope 
slides held at 84 K. After deposition, the film was ex- 
posed to the atmosphere for 10-20 minutes in order to 
allow a thin native oxide layer to form. Then a non- 
superconducting Al counterclcctrode was deposited from 
an Al 2024 alloy target, with the oxide serving as the tun- 
neling barrier. The low temperature parallel critical fields 
of the countcr-clcctrodcs were ^6 T, in good agreement 
with the expected i^c||■ The junction area was about 
1 mmxl mm, while the junction resistance ranged from 
15-100 kfi depending on exposure time and other factors. 
Only junctions with resistances much greater than that 
of the films were used, in order to be in the tunneling 
regime. 

II. VARIATIONAL BOGOLIUBOV-DE GENNES 
METHOD 

The combination of the Zeeman field and the disor- 
der potential ultimately leads to inhomogeneous, spin- 
dependent Hartree potentials. Therefore, we use a gener- 
alized Bogoliubov-de Gennes (BdG) method^ in which all 
2N BdG eigenvalues and eigenvectors are distinct (where 

is the number of sites). 

For convenience, we write the Hamiltonian in terms of 

an applied chemical potential /ir = /i — (where Vr is 
the quenched random potential) and field hr = h at every 
site: 

rr'(T r(T r 

(1) 

where Xra = fira — \ are densities with respect to half- 
filling and [/ < represents attraction. 

We decouple the Hubbard interaction in charge, spin, 
and pairing channels. It is difficult to justify a tradi- 
tional derivation of the self-consistent BdG equations 
with multiple-channel decoupling, because this appears 
to overcount the interaction term. We have performed 
a rigorous derivation based on the Trplnp variational 
formalism.^ In this approach, we postulate a trial Hamil- 
tonian Htriai, which dcfiucs a trial density matrix ptriai oc 
exp(— /3iftriai), and we then minimize the variational free 
energy ^ [given in Eq. (9)] with respect to the 3iV varia- 
tional parameters: the Hartree chemical potential /x^, 
Hartree field /i^, and self-consistent pairing field Af. 
This formalism has the practical advantage that can be 



used to assess the quality of the variational approxima- 
tion during the approach to self-consistency, and that it 
provides a rigorous upper bound to the true free energy. 
Our implementation is as follows: 

1. Make arbitrary initial guesses for the Hartree chem- 
ical potential /i^, Hartree field ^ and self- 
consistent pairing field Ar at every site r. These 
constitute a set of 3A'' real-valued variational pa- 
rameters. 

2. Find the effective chemical potential /Ir = Mr + 
and effective field /ir = /ir + at every site. These 
effective potentials include both the applied poten- 
tials and the Hartree potentials (resulting from the 
decoupling of the U term); they enter the mean- 
field Hamiltonian, 

H =^ Uy'cl^C^,^ -'^{Jir + hr(T)Xra- (2) 
rr'o" ra 

3. Construct the 2N x 2N Hamiltonian matrix 

Hra-.r'a': whcrc the indiccs (J, (t' =t:i distinguish 
between up-particle and down-hole sectors con- 
nected by matrix elements A: 

\ / acr' 

_,^flir + K _A A 

V Ar -Mr + KJ 

4. Diagonalize H to obtain eigenvalues Ea and eigen- 
vectors (/)arcr, where the eigenmode index a runs 
from 1 to 2N . (These eigenvectors are generaliza- 
tions of the and w^r vectors that appear in the 
original BdG formalism.) 

5. Find the symmetrized occupation numbers Ca = 
-itanhi;3E„. 

6. Compute the number densities x-^a (relative to half- 
filling) and the pairing density Fj. = i^^^c^^ at 
every site r: 

^rt = X] C,o.4>*oLvy4>a.v^, (4) 

a 

a^ri = - X] Ca</'ar;</'ar;, (5) 

OL 

= ^ Ca {4>l,yKvi + h.c) , (6) 

a 
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and thence the number density and magnetization 
on each site, 

Xr = i (Xr^ + Xrl) , (7) 
mr = i (Xrt -Xrl). (8) 

7. Compute the variational free energy 

O = ^ln(2coshi^i;«) 

a 

r 

+ J2'^i^rFr+IJ-^Xr + h^mr). (9) 
r 

8. According to the usual variational principle, we 
wish to minimize fl with respect to A, /i-^, and 

(to obtain a least upper bound to the true free 
energy). In practice this can be done by solving the 
stationarity condition V^^ = 0, i.e., finding a root 
of the SAT-dimensional equation 

f (X) = 0, (10) 

where X = {^^,h^,Ar} is the vector of varia- 
tional parameters and 

f = {^f + UXr, ftf - C/mr, Ar + U Fr} (11) 



is the residual vector (the "distance" from self- 
consistency). We use the standard Broyden 
method,^ which is a superlinearly convergent quasi- 
Newton method for multidimensional root-finding. 
The first iteration of the Broyden procedure is 
equivalent to fixed-point iteration of the self- 
consistency equations 

Af = -UFr, (12) 
Mf = -Uxr, (13) 
/if = +Umr. (14) 

We also inspect il to verify that the root of Eq. (10) 
corresponds to a minimum of Eq. (9) , and not to a 
maximum. We restart the Broyden method using 
Eq. (14) if a Broyden step results in a large increase 
in f2 (since quasi-Newton methods are prone to in- 
stability). 

After convergence we calculate further quantities, in- 
cluding the densities of states for up and down electrons 
(which are the main point of interest in this paper): 

a 

N,^{E) =Y,^iE + i?a)Cr40ari. (15) 
a 
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